Effective Medium Theory of Filamentous Triangular Lattice 
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We present a new effective medium theory that can include bending as well as stretching forces, 
and we use it to calculate mechanical response of a diluted filamentous triangular lattice. In this 
lattice, bonds are central-force springs, and there are bending forces between neighboring bonds on 
the same filament. We investigate the diluted lattice in which each bond is present with a probability 
p. We find a rigidity threshold independent of filament bending rigidity at pb — 0.56 and a crossover 
characterizing bending-, stretching-, and bend-stretch coupled elastic regimes controlled by the 
central- force isostatic point at pcF ~ 2/3 of the lattice when bending forces vanish. 

PACS numbers: 87.16. Ka, 61.43.-j, 62.20.de, 05.70.Jk 



I. INTRODUCTION 



Random elastic networks provide attractive and real- 
istic models for the mechanical properties of materials 
as diverse as randomly packed spheres [1~3], network 
glasses [4-8], and biopolymer gels [9-20]. In their sim- 
plest form, these networks consist of nodes connected 
by central-force (CF) springs to on average z neighbors. 
They become more rigid as z increases, and they typically 
exhibit a CF rigidity percolation transition [21-23] from 
floppy, disconnected clusters to a sample spanning-cluster 
endowed with nonvanishing shear and bulk moduli at a 
threshold very close to the Maxwell isostatic limit [24, 25] 
zcF = 2d, where d is the spatial dimension, at which the 
number of constraints imposed by the springs equals the 
number of degrees of freedom of isolated nodes. Gen- 
eralized versions of these networks, appropriate for the 
description of network glasses [4, 5] and biopolymer gels 
[13-15], include bending forces favoring a particular angle 
between bonds (springs) incident on a given node. For a 
given value of z, networks with bending forces are more 
rigid than their CF-only counterparts, and they exhibit 
a rigidity transition at ^ = Zb < zqf- 

Though numerical calculations, including the pebble 
game [23, 26], have provided much of our knowledge 
about the properties of random elastic networks, effec- 
tive medium theories (EMTs) [27-33] have provided sim- 
ple and qualitatively correct descriptions of CF networks. 
EMTs [34-36] and heuristic approaches [37] that describe 
both bending and stretching forces have only recently 
been developed. Here we present details of the deriva- 
tion of a bend-stretch EMT introduced in Ref. [35] and 
its application to a bond diluted triangular lattice, whose 
maximum coordination number is Zmax = 6. This lattice 
has bending and stretching forces modeled on those of 
biopolymer networks of filamentous semi-flexible poly- 
mers, characterized by one-dimensional stretching and 
bending moduli /i and k [13-15], respectively. Our EMT 
calculates the effective medium moduli u„, and 
function of /i and k and the probability p = z/Q that a 
bond is occupied. Both the EMT bulk and shear mod- 
ulus are proportional to When k = 0, our EMT 
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FIG. 1: (Color online) Phase diagram of the diluted filamen- 
tous triangular lattice showing the central-force and bending 
rigidity thresholds respectively at p = pcF and p = ph, the 
bending-dominated regime at small k in the vicinity of pb, 
the crossover bend-stretch regime near pcF, and the stretch- 
ing dominated regime at large k. 



reduces to that considered by others [29, 30] and suc- 
cessfully predicts a second-order CF rigidity threshold 
at zcF = 4 < Zmax [pcF = 2/3) with increas- 
ing linearly in p — pcF near pcF and approaching the 
undiluted triangular-lattice value of /i at p = 1. When 
bending forces are introduced, our EMT predicts a new 
second-order bending dominated rigidity threshold at 
P ^ Ph = 0.56 for all K > near which ^ '*(p ~ Ph) 
for K/{^a^) ^ 1 and /i^ ^ /J.(p — Pb) for K/{iJ,a^) ^ 1, 
where a is the lattice spacing. Near pcF, k is a relevant 
variable moving the system away from the CF rigidity 
critical point to a broad crossover regime [35, 37] in which 
Mm ^ K^^'^fJ.^^^ as shown in the phase diagram of Fig. 1. 
This crossover is analogous to that for the macroscopic 
conductivity in a resistor network in which bonds are 
occupied with resistors with conductance fT> with prob- 
ability p and with conductance cr< < o'> with probability 
l-p[38]. 

Though the model we study has both stretching and 
bending forces, it differs in important ways from previ- 
ously studied models for network glasses [4-8] and for 
filamentous gels [13-20]. The maximum coordination 
number for both of these systems is less than or equal 
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2d, and thus neither has a CF rigidity transition for 
p < 1 when there are no bending forces. As a result 
neither exhibits the bend-stretch crossover region near 
PCF that our model exhibits. Network glasses are well 
modeled by a randomly diluted four-fold coordinated di- 
amond lattice in which there is a bending-energy cost, 
characterized by a bending modulus k, if the angle be- 
tween any pair of bonds incident on a site deviates from 
the tetrahedral angle of 109.5 deg. The architecture of 
the undiluted diamond lattice (with z^ax = 4 < 2(i = 6) 
is such that its shear modulus vanishes linearly with k 
[7] and elastic response is nonaffine. When diluted, it 
exhibits a second-order rigidity transition from a state 
with bending-dominated nonaffine shear response to a 
state with no rigidity. As dilution decreases, rigidity is 
still controlled by k, but response becomes less nonaffine. 

Filamentous networks in two-dimensions are often de- 
scribed by the Mikado model [14-16] in which semi- 
flexible filaments of a given length L are deposited with 
random center-of-mass position and random orientation 
on a two-dimensional plane and the points where two fila- 
ments cross arc joined in frictionless crosslinks. As in our 
model, there is no energy cost for the relative rotation of 
two rods about a crosslink, but there is an energy cost 
for bending the rods at crosslinks. This model is charac- 
terized by the ratio rj = L/lc of the filament length L to 
the average mesh size, i.e., the average crosslink separa- 
tion Ic > a along a filament, where a is a minimum cutoff 
length. In the limit ?7 — )■ oo, all filaments traverse the 
sample, and the system has finite, K-indepcndent shear 
and bulk moduli: There is effectively a CF rigidity tran- 
sition at z = 4 when 77 is decreased from infinity. There 
is a transition at 77 « 5.9 to a rigid state with nonaffine 
response [14, 39], and there is a wide crossover region 
between rj — rjc and rj ^ 00 in which the shear modu- 
lus changes from being bend dominated, nonaffine, and 
nearly independent of /i at small 77 to being stretch domi- 
nated, nearly affine, and nearly independent of k at large 
77. Our EMT applied to the kagome lattice [40], whose 
maximum coordination number like that of the Mikado 
model is four, captures these crossovers. Interestingly, 
3d lattices composed of straight filaments with z,nax = 4 
exhibit similar behavior [41]. When filaments arc bent, 
however, elastic response in one case at least [20] is more 
like that of the diluted diamond lattice with the shear 
modulus vanishing with k even at at large L jl^ or z near 
4. 

External tensile stress (i.e., negative pressure) can 
cause a floppy lattice to become rigid [42]. Random in- 
ternal stresses can do so as well in a phenomenon called 
tensegrity [25]. Thus a lattice with internal stresses may 
have a lower rigidity threshold than the same lattice with 
out internal stresses [20] . Systems such as network glasses 
can exhibit two rigidity transitions [8, 43]: a second-order 
transition from a floppy to a rigid but unstressed state 
followed closely by a first order transition to a rigid but 
stressed state. These effects are beyond the scope of your 
EMT and will not be treated. 



The outline of our paper is as follows. Section II re- 
views properties of semi-flexible polymers and defines our 
model for the harmonic elasticity of crosslinked semi- 
flexible polymers on a triangular lattice; Sec. Ill sets 
up our effective medium theory; Sec. IV presents the re- 
sults of this theory; and Sec. V compares om- EMT with 
other versions of bend-stretch EMTs summarize our re- 
sults. There arc three appendices: App. B and App. C 
presents details of the derivation of crossover functions 
in the vicinity of pcF and pb, respectively; and App. D 
provides details on the estimation of pb using generalized 
Maxwell counting arguments. 



II. FILAMENTOUS POLYMERS ON A 
TRIANGULAR LATTICE 

The clastic energy of a single polymer is well described 
by the worm-like chain model [14, 44, 45] of an elastic 
rod 

where s is the arclength coordinate, L is the unstretched 
contour length of the polymer, dl{s)/ds and t(s) are, re- 
spectively, the extensional strain and unit tangent to the 
polymer at s, and ^ and k are, respectively, the one- 
dimensional stretching and bending stiffnesses. Three 
length scales can be identified in this clastic energy. The 
first is the con tour length of the polymers, L. The sec- 
ond, ^bcnd = \/ k/^, characterizes the relative strength of 
stretching and bending. For an elastic rod made of a ho- 
mogeneous material, /bend simply corresponds to the ra- 
dius of the rod. The third length, the persistence length, 
Ip = k/ (ksT), describes the scale of thermal undulations 
of the polymer at finite temperature. A fourth length, 
the mesh size Ic characterizing the connectivity of the 
network, can be identified for crosslinked polymer net- 
works. 

All these lengths participate in determining the elastic- 
ity of the polymer network above the rigidity threshold. 
Consider the force constants for harmonic distortions of 
an individual polymer segment of length Ic- At zero tem- 
perature, the force constant for longitudinal stretch of 
the segment is fcy^^ ~ m/^c, and that for transverse de- 
flection is fc_L ~ k//^ (ignoring numerical factors). Thus 
the ratio /bcnd/'c ~ V^-^/^W,^- characterizes the relative 
strength of bending and stretching of this particular seg- 
ment of length Ic at low temperature. At finite tem- 
perature, thermal undulations of the filaments induce an 
effective series spring constant fc|[ tj, ~ ksTlp/l'^ that re- 
sults from pulling out thermal undulations [13, 15, 17]. 
At any given temperature, the physical parallel inverse 
spring constant is simply k'^^ = fc|]~^ 4- k^li^ so that fcy ^b 
dominates at high temperature and ^ at low tempera- 
ture. The ratio Ic/lp ^ k±/k^l^ characterizes the relative 
ease of stretching vs. bending of this segment at high 



temperature. Finally, L/lc measures the connectivity of 
the network. 

Here we investigate the harmonic elasticity of a model 
filamentous network on a bond-diluted triangular lattice 
with lattice constant a which is the mesh length Ic- Poly- 
mers correspond to lines of connected, occupied colinear 
bonds and crosslinks to sites at which two polymers to 
crosslinks. We use a model elastic energy obtained by 
discretizing the contrinuum model of Eq. (1) on our tri- 
angular lattice. To simplify notation, we use ^/a = fi/lc 
to represent the full ki^\ rather than fcy ^^ only, and we set 
k± — n/a? = Thus the elastic energy for a defor- 

mation that maps the position of lattice sites ^ from 
to = Ti ~\- Ml can be written into the stretching part 
and the bending part 
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where u«/ = u^/ — u^, f«/ is the unit vector pointing 
from site I to nearest-neighbor (NN) site £' in the ref- 
erence state (i.e., the state with = for all £), gi (> 
is equal to one if the NN bond {£,£') is occupied and 
zero if the bond is empty. Owe" is the angle between 
bond {£, £') and {£' , £") in the deformed state as shown in 
Fig. 2(b), with r^, r^, and ri" along a line. The stretch 
spring constant /i/a is defined on each NN bond, and the 
bending constant n/a^ is defined for each next- nearest- 
neighbor (NNN) bond along a straight line (i.e., in the 
same filament). We have used the fact that the filaments 
are straight (6* = 0) in the reference state so f«' = f^'^". 
There is an important distinction, which will be of impor- 
tance in our development of an effective-medium theory, 
between the stretching energy Es and the bending energy 
Eh . The former is a sum of energies on individual bonds 
{£, £') connecting NN sites £ and £' , whereas the latter is 
a sum of energies on the pair of connected bonds {£, £') 
and {£',£") or, alternatively, on a phantom NNN bond 
connecting sites £ and f. This phantom bond is not 
an independent connection; it is only present and resists 
bending if both bonds {£,£') and {£',£") arc occupied. 
We will have more to say about this in Sec. V when we 
discuss the relation of our work to that of Refs. [341 and 



FIG. 2: (color online) (a) Filamentous triangular lattice with 
bonds randomly occupied with probability p. The unit vectors 
ej, e|, ef are marked by "1, 2, 3", the 3 stretch energy vectors 
are marked by the 3 red single lines, and the 3 bending 
energy vectors BJ^ are marked by the 3 green double lines. 
The purple dashed double line marks the bending vector Bl 
if the origin is marked by 0. (b) Illustration of the bending 
energy term 2b. The green curves correspond to the filament 
by requiring that it has to pass through lattice sites £,£',£". 
The ratio of the bending elastic modulus of the darker green 
segment and the lighter green segment is ki/k2 = 3. The 
bending energy term at site £' is proportional to 2(^ + ■^) 
determined by minimizing the elastic energy. 



III. EFFECTIVE MEDIUM THEORY 

We study the elasticity of our network using an 
effective-medium approximation [27, 28] that maps the 
random network to an effective-medium uniform lattice 
in which all bonds are present with respective stretching 
and bending constant fim and Km, which are determined 
self-consistently. This approximation has been shown 
to be a powerful tool for the calculation of properties 
of random systems, from the electronic structure of al- 
loys [27, 46] to the elasticity of random networks [29, 47]. 
It is straightforward to write the elastic energy AE of 
this uniform effective medium lattice in terms of the dy- 
namical matrix, Dq: 



AE = 



^ U_a • Da 
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where the Fourier transform of u is defined as 



Upon dilution, each of the bonds is present with a 
probability p, and the resulting lattice corresponds to a 
random network of semiflcxible filaments of finite ran- 
dom lengths L, whose average as a function of p is 
(L) = l/(l-p) [35]. 



and the dynamical matrix can be written as a sum of the 
outer product of the stretching and bending bond vectors 

3 K ^ 
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where 
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n,q 
6 

n,q 



= (l-e-'i-l^)ell 

- 2(l-cos(q.ell))ei^ 



(6) 



with e„ the unit vector along bond n defined in Fig. 2(a) 
and e;^ the unit vector perpendicular to elL The sum 
is over the three bonds in one unit cell for B'^ or the 
three pairs of bonds for B*". We used the periodicity of 
the lattice to obtain this translationally invariant form 
as will become clear when we discuss the perturbation 
calculation below. 

In the long- wavelength limit, the energy of our undi- 
luted harmonic triangular lattice reduces to the elastic 
energy of an isotropic 2d medium, 
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where u is the linearized symmetric Cauchy strain tensor 
with Cartesian components Uij, A and jl arc the Lame 
coefficients, X ~ fi = {\/3/A){fi/a), which depend only 
on fj, and not on k. The bending constant k only appears 
in the 0{q'^) terms in this elastic energy. 

The effective medium is characterized by a stretch 
spring constant that has the same value on every 
NN bond and a bending constant that has the same 
value on every NNN phantom bond along a line. The 
self-consistency conditions determining fi^ and k,„ are 
treated differently. In the dilution procedure, bonds are 
removed with probability 1 — p. In the process, a bond 
with spring constant Hm is replaced by one with spring 
constant 0, and the scattering potential relative to the 
effective-medium state arising from this replacement is 
proportional to /i— /i„i if the bond is occupied and — /i„i if 
it is vacant. Bending couples NNN sites along a line and 
must be treated differently. Consider two bonds along a 
line, {££') and {£'i"), sharing a common site £' and com- 
posed of rods with respective one-dimensional bending 
moduli Ki and K2 defined in Eq. (1). The phantom bend- 
ing bond connecting NNN sites is now a composite one 
[Fig. 2(b)] consisting of two NN bonds tied together at 
the intermediate site £' . A direct calculation of the min- 
imum bending energy [Eq. (1) with dl/ds = 0], for fixed 
bending angle 0, in this composite bond yields an effec- 
tive bending constant Koff = 2{n^^ + k^^)^^ that reduces 
to Ki if Ki = K.2 and if cither ki or K2 = 0. An example 
of the energy minimizing bending configuration is shown 
in Fig. 2(b). We derive the effective bending constant of 
this composite bond in App. A. 

Consider sites ^o, -^ii ^2, and £3 along a line. The re- 
moval of the central bond (Ixl-i) will modify the bending 
modulus of that bond and thus the bending energy of 
both NNN phantom bonds (IqI-i) and (£1^3) containing 
the NN bond {iii-i)- Thus, the bending modulus rela- 
tive to that of the effective medium of both these NNN 



phantom bonds is Kc('*s) — '^m, where 

,.1 IX- 
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where Ks = if the bond (^1^2) is absent and Ks = k if 
it is present. In the latter case Kc = 2KKm/{n + k,„) is a 
nonlinear function of k and Km- 

With these rules, we can write the perturbation to the 
dynamical matrix resulting from the replacement of one 
bond as 

Vq,q'(^s,Ks) = fl"^ (^is - Mm)Bl,qBi__q, 

+ a-3(K^(K,)-K„)Bj^qBj^„q, 

+ a-3(K,(K,)-K„0B^^qB^^_q,, (9) 

where we have chosen to replace the bond between the 
two sites ro = and ri = ae| for convenience. The 
vector B4 q = e~*'' ''i Bj q represents the bending of the 

pond pair connecting sites Tq = 0, ri = ae\, and 2ri, as 
shown in Fig. 2(a). The replaced stretching spring con- 
stant /is and bending modulus Ks follow the probability 
distribution 

P{Hs) = p6{^is - ^l) + {I ~ p)5{^is) 

P{ks) ^ p5{Ks-n) + {l-p)5{Ks). (10) 

In this study we are primarily concerned with the static 
properties of the effective medium, which is characterized 
by the zero-frequency phonon Green's function 



(11) 



When one bond of the effective medium is replaced, D — > 
D-|- V, and the Green's function is given by the following 
expansion 



Gl^' = (I-G.V)-lq,.Gq, 

— A^t'q,q'Gq + Gq • Tq q' • Gq' , 

where the T matrix is 

'^q.q' — ^q,q' + T7 5^^q,qr^qi'^qi,q' 



(12) 



qi 



~'~7V2 ^q,qi'^qi'^qi,q2'^q2'^q2 q' 

qi,q2 

+ •••. (13) 

In EMT, the effective- medium elastic parameters //^ and 
^771 are determined by requiring that the disorder average 
of G^ be equal to G of the effective-medium uniform 
lattice with stretching and bending constants and Km- 
This is equivalent to requiring that the disorder average, 

pT|^^=^,,^=« + (1 - p)T|^^=o,«,=o = 0, (14) 
of the T-matrix vanish. 
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In traditional EMT with no bending modulus, the only 
elastic parameter is the stretching spring constant 
and one can simplify the T matrix by making use of the 
equation 

qi 

where / is a scalar function. However, in this case with 
bending energy terms, which couple the replaced bond 
with two neighboring bonds, the product Vq qj-Gqj-Vq^ .q' 



contains cross terms such as B? „ -G^.-B^^ B^J „, 

which do not reduce to Vq.q//. 

Here we develop a new variation of the EMT to solve 
this self-consistency equation (14). In this method, we 
introduce an effective medium coupling between neigh- 
boring bond pairs along filaments that share a bond, e.g., 
BJB4. with a coupling constant Am, to account for the 
cross terms discussed above. For convenience we express 
the self-consistency equation in the space defined by the 
basis {|Bf), [B^), |B|)}, in which the perturbation is 



(P-s 








^)/a 



(16) 



where the tilde marks the addition of the Am coupling. 
Cross terms such as BfBj vanish because they violate 
the space inversion symmetry, q — )■ — q. 

On the other hand, the Am coupling terms also need 
to be added to the dynamical matrix 

-Dq(/^?n7 ^rn) -Dq(/im, l^ni) 
3 

+a-3Am^2cos(q-ell)B^,qB:;__q, (17) 

n=l 

and the Green's function is modified to G = — D^^. 
By writing the T matrix as 

t=(V-i-G)^\ (18) 

the self-consistency equation takes the form 

p[V-\fi,n) - G]-i + (1 -p)[V-i(0,0) - G]-i = 0(19) 

Multiplying by (V~^(/i, /c) — G) on the left and by 
(V~^(0,0) — G) on the right we arrive at 

pV-\0, 0) + (1 - p)V-i(/i, k) - G = 0. (20) 

This equation can be solved in the space of 
{|B^), |B5), |B|)}. To do this we define the notation 

H(6m, Im) = - — G(^m, Km, Am)- (21) 

a 

From the definition of G it is straightforward to see that 
H only depends on the ratios bm = Km/fJ-m and Im = 
Am/Mm- We then project H to this basis 

/ Hi \ 
H(&m,Z™) = i/2 i/3 , (22) 
V H3 H2 J 



with 

Hi^{Bl\li\Bl)^ f ^B-J.q.Hq.B^ (23) 

where we have defined the inner product as an integral 
over the first Brillouin zone and the factor v is the volume 
of each unit cell in real space. Similarly, H2 = (BJ|H|B5) 
and H3 = (B5IHIB4). As we mentioned earlier, cou- 
plings between B** and B^ terms vanish in the dynamical 
matrix D, because D is even in momentum q. As a re- 
suh, the projections (Bf |H|B5) = and (Bf |H|B^) = 0. 
This is consistent with the vanishing cross term between 
B^ and B^ in V~^. Thus the self-consistency equa- 
tion (20) closes, providing 3 independent equations for 
the 3 unknowns (pLrm ^m, Am) whose solution we will dis- 
cuss in the next section. 



IV. EMT RESULTS 

Numerical solutions to Eq. (20) for any value oi k/h 
are easily calculated, and the results for the effective 
medium elastic parameters, and the results are plotted 
in Figs. 3, 4. There are several properties of these plots 
that are worthy of note: 

1. /im. vanishes at the CF Maxwell rigidity threshold 
PCF = 2/3 when k = and at p = = 0.56 for 
all K > 0. Simulations of the same model yield 
that PCF ^ 0.659 and pb ^ 0.445 [36]. (In App. D, 
we present a variation of the Maxwell floppy mode 
count to estimate the rigidity threshold in presence 
of filament bending stiffness and obtain pb — 0.448 
in good agreement with simulation results.) 

2. Increasing k increases fim for all p except at p = pb- 

3. For small K/{^a'^), there is an interesting and 
nontrivial crossover near pcF whereas for large 
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K/{^a?), memory of the CF threshold is effectively 
lost and rapidly reaches value near its satura- 
tion value /i for p > p^. 

4. Km vanishes as p — > pb and rises smoothly to its sat- 
uration value K without any evidence of crossover 
behavior near pcF • 

5. Xrn vanishes at pb and in the undiluted lattice 
{p ~ 1), which it must by construction. It exhibits 
crossover behavior near p = pcF for small K/(/ia^). 

In what follows, we discuss analytic expressions derived 
from Eq. (20) that captures these behaviors. 
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FIG. 4: (color online) The EMT solution for (a) and A„ 
(b). Parameters and color code are the same as in Fig. 3. 



of K as 



C25p 



a? K/{^a^) + Ci ' 



(24) 



0.6 0.7 
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(b) 

FIG. 3: (color online) (a) Semi-log plot for the EMT 
solution 

fJ/rn SI'S db function of p for jj. = 1 and k = 
1,10"^10~^1Q~^10"*,10"^10"'' from top to bottom, as 
indicated in the legend, (b) Linear plot of /im as a function of 
p, with parameters and color code the same as in (a). The red 
solid line indicates fim for the case of a central force triangular 
lattice (k = 0). 



A. Behavior near p = pb 

We begin with behavior near pb- As detailed in App. C, 
firn approaches zero linearly in Sp = p — pi, for all values 



where ci and C2 are constants determined by the ar- 
chitecture of the lattice and are independent of p or 
K/(//a^), and both Km and Xm are linearly proportional 
to fj-m for small Sp with corrections of order Sp^. The 
linear dependence on Sp of these quantities is what one 
would expect from and EMT. The prefactor, shows im- 
portant dependence on k. When K/{^a?) <C 1, /im is 
linearly proportional to k and independent of whereas 
for K/{^a^) 3> 1, /im is linearly proportional to fi and 
independent of k. Thus the system behaves as if the it 
consists of a "bending spring" with spring constant k/u^ 
in series with a stretching spring with spring constant 
fi/a. This result explains the increase in slope of at 
Pb with increasing k and the saturation of for p near 
Pb for large K/{^a?'). In Fig. 5 we show a comparison 
between the asymptotic solution (24) and the numerical 
solutions of the EMT equations, which agree very well. 
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FIG. 5: (Color online) Asymptotic solution (dashed lines) and 
numerical solutions (data points) of fim near ph- Parameters 
and color code are the same as in Fig. 3. 



B. Behavior near p — pcF 

We show in App. B that in the vicinity of the CF rigid- 
ity threshold at small k, the effective medium spring con- 
stants for given Ap = p — pcF adopt the scaling forms 



Mm = M 



- A^|Ap^^92,± 



K 



(25a) 
(25b) 



,a2/i|Ap|'^/' 
A„ = M|Apr^.3,.(^^^), (25c) 



where 





= 2, 


(26a) 


h 


= 1, 


(26b) 




= 2, 


(26c) 


h 


= 3, 


(26d) 



and 




92,±{x) ~ -X, 



(27a) 
(27b) 



g^^^[x) ^ ^\±l + \ll-^x\ x\ (27c) 



where A ~ 2.413 and B ~ 1.520 arc constants evalu- 
ated in App. B. These scaling relations are analogous to 
that found in random resistor networks with two different 
types of resistors [38] , and central force spring networks 
with strong and weak springs [37]. 

Plots of the numerical results in these scaling forms are 
shown in Fig. 6 along with the above analytical result. 
These scaling relations can also be expressed in terms of 
the mean filament length (L) = 1/(1— p) thus (L)cF = 3 



and|p— pcfI (X — (L)cf|. Thus, the scaling exponents 
for functions expressed in terms of (L) are the same as 
those expressed in terms of p. 

Macroscopic elasticity of the effective medium can be 
described by the Lame coefficients, which are both lin- 
early related to /im, as we discussed in Eq. (7). Thus, the 
asymptotic forms of /im can characterize different elastic 
regimes of the diluted lattice: 



MiApi; 



± 1 



'1- 



9 a2^|Ap|'? 



(28) 



3/x|Ap| 

A K. 

V 3a^ |Ap| 



»1, 



- < 1 and Ap > 0, 
< 1 and Ap < 0. 



These crossover regimes correspond exactly to those 
found in Rcf. [37] using known behavior of the density 
of states and mode structure of systems near the CF iso- 
static limit and general scaling arguments. 

Thus, it is straightforward to see that above the 
crossover regime, the macroscopic elastic moduli are pro- 
portional to /i, indicating that the elasticity of the net- 
work is dominated by stretching energies, whereas below 
the crossover regime, the moduli are proportional to k, 
indicating that the network elasticity is dominated by 
bending energies, in agreement with the fact that the 
system is below the central force isostatic point. Interest- 
ingly, in the crossover regime near pcF, the macroscopic 
elastic moduli is proportional to p}/'^K^/'^ , corresponding 
to strong coupling between both stretching and bending 
modes. 

These observations indicate a phase diagram as shown 
in Fig. 1. There is a rigidity transition between the state 
with no rigidity (vanishing elastic moduli) and rigid state 
at pb- Above this rigidity transition and at small filament 
bending stiffness k, the elasticity of the diluted filamen- 
tous lattice is controlled by a crossover at the central force 
isostatic point pcF- For large bending stiffness (k ^ 
the system directly enters the stretching dominated affine 
elastic regime as p is increased. 



V. DISCUSSION 

Two other approaches produce results similar to ours, 
and we briefly review compare them to ours. 

In references [34] and [36] , Das et al. present an alter- 
native EMT for systems with both stretching and bend- 
ing forces. The first of these reference fails to identify 
either the new rigidity threshold at pi, and the wide 
bend-stretch crossover regime near pcF, but the sec- 
ond does. Stretching forces are easily described by CF 
springs, which reside on bonds, each of which can have 
a distinct spring constant, whereas bending forces effec- 
tively reside on phantom NNN bonds. Because remov- 
ing one NN bond from a pair defining a phantom NNN 
bending bond effectively removes that phantom bond. 
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10' 




(P-Pcf)* 

(a) 




0.001 0.1 10 1000 



K 

(P-Pcf)* 

(b) 




10-* 0.01 1 100 10^ 



K 

(P-Pcvf 

(c) 

FIG. 6: (color online) (a), (b), and (c) show rescaled plots 
of the EMT solutions /im, Km and Am using the scaling 
forms (25a) for ^ = 1 and k. = 10"^ 10"^ 10"*, 10"^ 10"^ 
with color code the same as in Fig. 3, and exponents taking 
the value as in Eq. (26a). The thin black lines represent the 
asymptotic forms (27a) for small k as solved in App. B. The 
brown dash-dotted lines, the thick blue solid lines, and the 
purple dashed lines plot the functional form of jim obtained 
in the crossover, the stretching-dominated, and the bending- 
dominated regimes of Eq. (28) , respectively. 



bending and stretching are not independent in the di- 
luted lattice. This presents real challenges for the devel- 
opment of a consistent bend-stretch EMT. Our approach 
to this problem appeals to the underlying polymer na- 
ture of our model in which constituent polymers are en- 
dowed with local stretching and bending moduli ^ and k. 
We can modify these moduli along any bond. Different 
stretch moduli lead to independent effective CF stretch 
force constants fc|| for each bond. Modification of the 
bending modulus k on a given NN bond, however, mod- 
ifies the bend force constant for both phantom NNN 
bonds that that NN bond partially defines in the manner 
described above. With this approach, we develop a con- 
sistent EMT that includes the statistical correlation be- 
tween bend and stretch. Das et al. ignore this correlation 
and assume that a stretch spring on a given NN bond can 
be removed without affecting the bending energy on the 
phantom NNN bonds that include that NN bond and 
that bending springs on the phantom NNN bonds can 
be removed without affecting the stretch springs on the 
two bonds that define the NNN bond. In other words, 
the phantom bond is effectively elevated to a real bond 
with existence independent of the underlying NN bonds. 
In general, the NNN bonds can be present with with an 
arbitrary probability q and absent with probability 1 — q. 
To provide an approximate description of the constraint 
that the phantom bond does not exist unless both of 
the NN bonds defining it are present. Das et al. assign 
a probability q = [p is the probability that a NN 
bond is occupied) to the occupancy of a NNN bending 
bond, but continue to treat the NN and NNN bonds 
as statistically independent. Again the result is a set of 
closed self-consistent equations for fim and Km- Both ap- 
proaches yield pcF = 2/3 in good agreement with simula- 
tions [35], which yield pc¥ = 0.651. Our approach yields 
a value for p^ (0.56) that is well above that (0.445) ob- 
served in simulations [35] whereas that of Rcf. [36] yields 
a value (pb = 0.4) in good agreement with simulations. 
Both approaches yield a nontrivial bend-strctch crossover 
in the vicinity of pcF in qualitative agreement with simu- 
lations. Our approach gives values of the shear modulus 
in the crossover region of the same order of magnitude as 
seen in simulations whereas that of Ref. [36] finds values 
1.5 to 2 orders of magnitude higher at small values of n. 

In Ref. [37], Wyart et al. consider random off-lattice 
elastic networks derived from two-dimensional packings 
of spheres [48] with a cordination number above the 
Maxwell CF isostatic limit of z 4 in which CF springs 
are assigned to each sphere-sphere contact. They use 
numerical simulations to study the nonlinear relation be- 
tween shear stress a and shear strain 7 as springs are 
cut, thereby reducing z, and they find a scaling relation 
a = 7[(5z[/(7/[(5z[), where (5z = z — 4, f{x) — >■ const, for 
X — > 0+, f{x) — > for .T — > 0~, and f{x) ~ x for x — > 00. 
This scaling form predicts (7^7 for Jz3>7>0,(t = 
for 5z <C —7 < 0, and cr '--^ 7^ for 7 ^ 5z. Reference [37] 
then provides a theoretical justification for this behavior 
based upon the existence of a plateau in the density of 
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states [49] above to* ~ 5z and reasonable assumptions 
about statistical independence of eigenvectors associated 
with different normal modes in the isostatic network [50] 
and about the nature of nonaffine response of nearly iso- 
static systems. Finally, they extend this line of reasoning 
to nearly isostatic systems with extra weak bonds and 
find three regimes of elastic response that arc identical 
to those we identify in Eqs. (25a) to (27a) if the weak 
bonds are of a bending type. 

To summarize we developed a new effective-medium 
theory that can include bending energy of filaments, 
and we used it to study the development of rigidity 
of a randomly diluted triangular lattice with central 
force springs on occupied bonds and bending forces be- 
tween occupied bond pairs along a straight line. We ob- 
tained a rigidity threshold for positive bending stiffness 
and a crossover, controlled by the isostatic point of the 
central force triangular lattice, characterizing bending- 
dominated, stretching-dominated, and stretch-bend cou- 
pled clastic regimes. 

Acknowledgments — We are grateful to C. P. Broedersz 
and F.C. MacKintosh for many stimulating and helpful 
discussions. This work was supported in part by NSF- 
DMR-0804900. 



Appendix A: Bending of a composite rod 

In this Appendix we solve for the effective bending 
rigidity of a composite rod. We solve this problem using 
a simple setup shown in Fig. 7. We assume that there is 
a bond of bending rigidity ki connecting points ££' , and a 
bond of bending rigidity K2 connecting points £'£" , with 
points £,£',£" hold in positions 



{-1,/io}, {0,0}, {i,M, 



(Al) 



respectively in space, without losing generality. We then 
calculate the shape of the bonds which minimize the elas- 
tic energy given these fixed end points, treating the bonds 
as elastic rods. In addition we assume that the stretching 
stiffness is much greater than the bending stiffness so we 
can ignore the length change of the bonds. The elastic 
energy of a bond is then simply the bending energy of a 
rod with bending stiffness k 



E 



dx 



d^h{x) 
dx^ 



(A2) 



where we have assumed small deformations so the elastic 
energy takes the form of an expansion in small h{x). The 
variational condition that a configuration minimizes this 
elastic energy is 



d^h{x) 



= 0. 



(A3) 



Therefore, we can assume the shape of the two bonds 

h(x) ^ Ax + -Bix^ + -Cix^ for a; < 
2 6 

h(x) ^ Ax + ^B2X^ + \C2X^ for x > (A4) 
2 6 

where the first line corresponds to the bond between ££', 
and the second line correspond to the bond between £'£". 
This form already satisfies the condition that the two 
bonds are connected and smooth at their conjunction 
X = 0. We have two boundary conditions at a; = — 1 and 
X = 1 



ho 
ho 



-A 



A 



2 6 

— Bo -\ — Co . 
2 6 



(A5) 



Under these boundary conditions we minimize the total 
elastic energy of the two bonds 



E^^ dx 
2 J-i 



'd'^hix)' 


2 

H 




'd'^h{x)' 


dx^ 






dx'^ 



1^2 



i?2 + B2C2 



CI 



2 V 3 
This calculation leads to the minimum clastic energy 

6KlK2/lg 



(A6) 



Em. = 



Kl + K2 



(A7) 



In the above calculation, if we assumed that both of 
the two bonds are of bending stiffness K^g, we get 



E„ 



(A8) 



Therefore we arrive at the effective bending stiffness of 
the composite rod 



2kiK2 



Kl + K2 

which reduces to ki if ki = K2 as we expected. 



(A9) 



Appendix B: Asymptotic solution of the CPA self 
consistency equation near pcF 

In this Appendix wc solve the EMT self-consistency 
equation [Eq. (20)] asymptotically near pcF at small k. 
In the zeroth order, k = 0, and the problem reduces to 
that of a central force rigidity percolation. We expand 
the equation around this point and solve for the effective- 
medium spring constants asymptotically in small k. 

To obtain the asymptotic solution, we rewrite the EMT 
self-consistency equation [Eq. (20)] as 

pYiii, k) + (1 - p)V(0, 0) - V(0, 0)GV(/i, k) = 0,(B1) 
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r 



FIG. 7: (Color online) Bending of the composite rod. The 3 
red disks mark the fixed end points of the two bonds. The 
blue solid line represent the shape of the two bonds solved 
for the case of ki — 10k2, and the red dashed line represent 
the shape of the two bonds for the case of ki — K2 which is 
symmetric under x — > — a-. 



to avoid the apparent singularity in V ^ at the point 



K = 0. 

It is straightforward to see from the symmetry of 
and G as in Eqs. (16) and (22) that the matrix equation 
[Eq. (20)] contains 3 independent equations. The 11- 
element gives 



p ~ Hi{brn,lm) 
1 — Hi{bm,lm) 



(B2) 



To zeroth order, k = 0, we obtain iJi(0,0) = pcF = 
2d/z ~ 2/3 from the fact that all bonds are equivalent, 
and Eq. (B2) reduces to the EMT equation for central- 
force rigidity percolation. 



The 22- and 23-clements of the matrix equation (20) read 



2(i + 5^ 



2'i + 5^ 



bmM?. — Im.H? 



(52 



(lmH2 + bmH'i) + 2b„il„iH2 



lfjH2 + 26„/™i/3 = 0, 
(6^-1-4)^3 = 0, 



(B3) 
(B4) 



where b = K/{fima^), H2{brn,lm) and H^{b,n,lm) are the 22 and 23 elements of H{b,n,lm) respectively. Thus we 
have 3 imknowns {/im,6m, /„i} (or equivalently, {/im, k™, A„i}) and 3 equations (B2), (B3), and (B4). As we already 



discussed, at k = 0, the zeroth order solution is {/^ 



(0) 



PCF 1,(0) 



1-PCF ' 



0,1 



(0) 



0}, where PCF = -ffi(0,0). Notice 



that the equations up to this point are all exact within the EMT. 

As K becomes positive, increases, bm and Im become nonzero, and the rigidity threshold jumps to a lower value 
Pb as shown in Fig. 3(b). For small k, we have k/I^o^) ^ 1, and can assume that bm,lm ^ 1 (which we will verify 
later), and we find that to the leading order the 3 Eqs. (B2), (B3), and (B4) become 



P-PCF - i?l,l(0, 0)K„/(/Xrna^) 



1 -PCF 



k{2p - 1), 
Ki^3(0,0)- 



K 1 



(B5a) 
(B5b) 
(B5c) 



where i?i,i(0,0) = dHi/dh 



m b„,=0,Z„=0 



MmO" P 

—2.413 and i?3(0, 0) = 1.520. From these relations, we find that at p = pcF 
K^/^, (B6a) 



Am 



.3/2 



(B6b) 
(B6c) 



indicating that ^ Km ^ Am and thus b„i,lm ^ 1 as we assumed. Using these relations, together with the fact 
that as K — >■ 0, fim P ^ Pcf, we arrive at the scaling relations 



K 



K 



(B7a) 
(B7b) 
(B7c) 



where Ap = p — pcF- These scaling forms agree well with our numerical solutions (See Fig. 6). 



11 



From Eq. (B5), we can further solve for the form of these scaling functions 



9i..± {^) - \ \±l + \jl- ^xj , (B8a) 

g2,±{x) ~ ^x, (B8b) 

-1 



g3,±(x) :^ ^ ±l + ^l-i^:,l (B8c) 



with ^ = -Fi,i(0,0) and e = i?3(0,0). These scaling 
forms with Hi i and obtained numerically are plotted 
together with the numerical solutions in Fig. 6 in nice 
agreement. 

These scaling functions can be expanded in different 
limits. In particular, for the effective medium stretching 
stiffness ^, which is closely related to the macroscopic 
elastic moduli, we have 

{\fAy/x ifx>l, 
3 if X < 1 and Ap > 0. (B9) 

\Ax if x < 1 and Ap < 0, 

indicating different elastic regimes as we discussed in 
Sec.III. 



Appendix C: Asymptotic solution of the CPA self 
consistency equation near pb 

Equations (B2,B3,B4) can also be used to solve for the 
asymptotics near the rigidity threshold ph- In particu- 
lar, because lm,bm converge to constants that are much 
smaller than unity and independent of k near p^, the 
asymptotic solution near pb in this section are not lim- 
ited to small K. 

Firstly we solve for the value of the rigidity threshold 
Pb for the case of k > using these EMT equations. At 
Pb, we have fim = 0, k,,,. = 0, Am = and as a result 
6 — > oo, and we solve for Ph,bm,lm- So the equations 
determine pb,bm = bb,lm = h are 

Pb - Hi{bb,lb) = 0, 

26t(pb - ^) + i-bl + ll)H2{bb,lb) = 0, 

-lb + {-bl + ll)H3{bbJb) = 0, (CI) 

where bb and lb are the value of and Im at pb. This 
set of equations is independent of k. Numerical solution 
to these equations are given by 



Ph 



lb 



0.5584, 

0.06355, 

0.004235, 



(C2) 



which agrees with the results we obtained by solving the 
EMT equations numerically. 



Secondly we solve for the asymptotic behaviors near 
Pb. To achieve this, we suppose p = pb -I- 5p, and to first 
order we have 



Mm = + (5yU„ 

bm = bb + Sbr, 

Im lb ^ ^Im 



(C3) 



We put these expansions back into Eqs. (B2,B3,B4) we 
get the first order perturbation equations 



^/i,n = A^ '"^""^r"p: ' (C4) 
(2pb - 1 - 2bbH2.o)5b„, = -2bbSp + 2pb6g%^,(C5) 

(C6) 



SL 



-2bbH-- 



In deriving these equations we used the fact that lb <C 
bb I and ^ = ^ ^ 1 near pb. Thus we arrive at the 
asymptotic solution of the effective medium stretching 
stiffness 



C2<5p 



(C7) 



where 



Cl 



C2 = 



Ai-2bbH3,oA2 2phbl 
2pb - 1 - 2bbH2.o 1 - Pb 

1 Ai-2bbH3nA2 



1 -Pb 



l + 26fa 



2pb - 1 - 2bbH2.o 



(C8) 



are constants determined by the architecture of the lat- 
tice and are independent of p or k/Iho^). In the case of 
triangular lattice we have ci = 0.1018 and C2 = 5.132. 



Appendix D: Estimates of the rigidity threshold 
using Maxwell counting of floppy modes 

In this appendix, we use Maxwell counting [5, 6, 24] 
of the number of fioppy modes to estimate the rigidity 
threshold of a filamentous triangular lattice. Bending 
forces cause bending excitations on a single filament to 
have nonzero frequency. Thus in the counting of zero 
modes, we can treat each filament as a rigid rod with 
two translational and one rotational degree of freedom. 
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Thus the number of floppy modes in the network is given 

by 



(Dl) 



where Nf is the total number of filaments and N^. is the 
total number of constraints in the network. Upon dilu- 
tion, Nf increases and Nc decreases, and thus iVo be- 
comes positive at the threshold occupancy probability 
Ph- 

To calculate the number of filaments Nf, we use the 
fact that the average length of the filaments is [35] 



1 



1-p 



(D2) 



Consider a straight line in the lattice (which would be 
one whole filament if p = 1) with Nj, bonds, the mean 
number of filaments in this line Nf satisfies the equation 



Nf{L)^pNb. 



(D3) 



Thus we get the mean number of filaments per lattice site 



^3p{l-p). 



(D4) 



The number of constraints at a crosslink where 2 fil- 
aments meet is 2, because the crosslink puts the con- 
straint on the 2 filaments that the linked sites on these 

2 filaments have to move together in both x and y di- 
rections. The number of constraints at a crosslink where 

3 filaments meet is 4 correspondingly. Thus for occu- 
pancy probability p of the bonds, the mean number of 
constraints per lattice site is given by 

Uc = 4p*^ + 4 X 6p^{l - p) 
+(4 X 12 -f 2 X 3)/(l 
+(4 X 8 + 2 X 12)p^{l-pf 
+2xl2p\l-p)'^, (D5) 

with the terms illustrated in Fig. 8. 

From the equation of fioppy modes per lattice site 



3nf — He ^ 
we obtain the rigidity threshold 
Pb = 0.4475, 



(D6) 
(D7) 



which is very close to the value obtained in numerical 
simulation of filamentous triangular network [35] . 

This estimation of the rigidity threshold can also be 
applied to other cases. For example, we consider the 
two dimensional kagome lattice, which has z = 4 when 
undiluted. Thus the p ^ 1 point is the CF isostatic 
point. In the presence of bending forces, we follow a 
similar calculation and get 



nf 



2p{l 
2/4 



2 X 4p3(l 




A 




is) 



FIG. 8: Illustration of the configurations of bonds connecting 
to a site for given probability p, corresponding to terms in 
Eq. (D5) (only one diagram is shown for each type of config- 
urations with the same probability and the same number of 
constranits). From (a) through (g) the probability and the 
number of constraints are respectively (p*,4), (6p^(l — p),4), 
(12p^(l-p)^4), (3/(l-p)^2), (8p«(l-p)^4), (12p=*(l- 
pY,2), (12p^(l -p)'',2). The rest of the possible configura- 
tions, such as only two bonds along a line or only one bond 
connected to the site does not contribute constraints. 



Thus the rigidity threshold with bending forces is at 

Pb = 0.6074. (D9) 

Another example we consider is the three dimensional 
filamentous face centered cubic (FCC) lattice. The undi- 
luted FCC lattice has z = 12 and is well above isostatic- 
ity, similar to the case of triangular lattice, n/ and tlc 
for the FCC lattice is given by 



Uf 



-p) + 2 X 4p2(i-p)2. (D8) 



6p(l -p), 

3[5p^2 + 5x 12p"(l-p) 

-f (5 X 60 + 4 X 6)/°(l -p)2 

+ (5 X 160 + 4 X m)p^{l-pf 

+ (5 X 240 + 4 X 240 + 3 X 15)/(1 - pf 

+ (5 X 192 + 4 X 480 + 3 X 12G)/(1 - pf 

+ (5 X 64 + 4 X 480 + 3 X 360 + 2 X 20)/(l - pf 

+ (4 X 192 + 3 X 480 + 2 X 120)/(1 - pf 

+ (3 X 240 + 2 X 240 + 1 X 15)/(1 - pf 

+ (2 X 160 + 1 X m)p^{l-pf 

+60p2(l -p)i"]. (DIO) 

In three dimensions the equation for vanishing floppy 
modes is 

5n/ - = (Dll) 
and we get the rigidity threshold for FCC lattice at 

Pb ^ 0.3062. (D12) 



13 



[1] A. J. Liu and S. R. Nagel, Nature 396, 21 (1998). 

[2] M. Wyart, Ann. Phys. Fr 30, 1 (2005). 

[3] A. J. Liu and S. R. Nagel, in Annual Review of Condensed 

Matter Physics (2010), vol. 1, pp. 347-369. 
[4] J. C. Phillips, J. Non-Cryst. Solids 43, 37 (1981). 
[5] M. Thorpe, J. Non-Cryst. Solids 57, 355 (1983). 
[6] J. C. Phillips and M. F. Thorpe, Solid State Commun. 

53, 699 (1985). 
[7] H. He and M. F. Thorpe, Phys. Rev. Lett. 54, 2107 

(1985). 

[8] M. Thorpe, D. Jacobs, M. Chubynsky, and J. Phillips, J. 

Non-Cryst. Sohds 266, 859 (2000). 
[9] E. Elson, Annu. Rev. Biophys. Chem. 17, 397 (1988). 
[10] K. Kasza, A. Rowat, J. Liu, T. Angelini, C. Brangwynne, 

G. Koenderink, and D. Weitz, Curr. Opinion Cell Biol- 
ogy 19, 101 (2007). 
[11] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, 

and P. Walter, Molecular Biology of the Cell (Garland, 

New York, 2008), 4th ed. 
[12] P. Janmey, S. Hvidt, J. Lamb, and T. Stossel, Nature 

345, 89 (1997). 
[13] F. C. MacKintosh, J. Kas, and P. A. Janmey, Phys. Rev. 

Lett. 75, 4425 (1995). 
[14] D. A. Head, A. J. Levine, and F. C. MacKintosh, Phys. 

Rev. Lett. 91, 108102 (2003). 
[15] J. Wilhelm and E. Frey, Phys. Rev. Lett. 91, 108103 

(2003). 

[16] D. A. Head, A. J. Levine, and F. C. MacKintosh, Phys. 

Rev. E 72, 061914 (2005). 
[17] C. Storm, J. Pastore, F. MacKintosh, T. Lubensky, and 

P. Janmey, Nature 435, 191 (2005). 
[18] P. R. Onck, T. Koeman, T. van Dillen, and E. van der 

Giessen, Phys. Rev. Lett. 95 178102(2005). 
[19] E. M. Huisman, T. van Dillcn, P. R. Onck, and E. Van der 

Giessen, Phys. Rev. Lett. 99, 208103 (2007). 
[20] E. M. Huisman and T. C. Lubensky, Phys. Rev.Lett 106, 

088301 (2011). 

[21] S. Feng and P. N. Sen, Phys. Rev. Lett. 52, 216 (1984). 
[22] S. Feng, P. N. Sen, B. I. Halperin, and C. J. Lobb, Phys. 

Rev. B 30, 5386 (1984). 
[23] D. J. Jacobs and M. F. Thorpe, Phys. Rev. Lett. 75, 

4051 (1995). 
[24] J. C. MaxweU, Philos. Mag. 27, 294 (1864). 
[25] C. Calladine, Int. J. Solids Struct. 14, 161 (1978). 
[26] D. J. Jacobs and M. F. Thorpe, Phys. Rev. E 53, 3682 

(1996). 



[27] P. Soven, Phys. Rev. 178, 1136 (1969). 

[28] R. J. Elliott, J. A. Krumhansl, and P. L. Leath, Rev. 

Mod. Phys. 46, 465 (1974). 
[29] S. Feng, M. F. Thorpe, and E. Garboczi, Phys. Rev. B 

31, 276 (1985). 

[30] E. J. Garboczi and M. F. Thorpe, Phys. Rev. B 31, 7276 
(1985). 

[31] L. M. Schwartz, S. Feng, M. F. Thorpe, and P. N. Sen, 

Phys. Rev. B 32, 4607 (1985). 
[32] C. Heussinger and E. Frey, Phys. Rev. Lett. 97, 105501 

(2006). 

[33] C. Heussinger, B. Schaefer, and E. Frey, Phys. Rev. E 

76, 031906 (2007). 
[34] M. Das, F. C. MacKintosh, and A. J. Levine, Phys. Rev. 

Lett. 99, 038101 (2007). 
[35] C. P. Broedersz, F. MacKintosh, X. Mao, T. Lubensky, 

and F. MacKintosh, arXiv:1011.6535vl [cond- mat. soft] 

(2010) . 

[36] M. Das, D. Quint, and J. Schwarz, arXiv:1106.3004 

(2011) . 

[37] M. Wyart, H. Liang, A. Kabla, and L. Mahadevan, Phys. 

Rev. Lett. 101, 215501 (2008). 
[38] J. P. Straley, J. of Phys. C Solid State 9, 783 (1976). 
[39] M. Latva-Kokko and J. Timonen, Phys. Rev. E 64, 

066117 (2001). 

[40] X. Mao, O. StenuU, and T. C. Lubensky, Unpublished 
(2011). 

[41] O. StenuU and T. C. Lubensky, arXiv:1108.4328vl [cond- 

mat.soft] (2011). 
[42] S. Alexander, Phys. Rep. 296, 65 (1998). 
[43] P. Boolchand, G. Lucovsky, J. C. Phillips, and M. F. 

Thorpe, Philos. Mag. 85, 3823 (2005). 
[44] O. Kratky and G. Porod, Rec. Trav. Chim. Pays-Bas. 68 

(1949). 

[45] J. Marko and E. Siggia, Macromolecules 28, 8759 (1995). 
[46] S. Kirkpatrick, B. Velicky, and H. Ehrenreich, Phys. Rev. 

B 1, 3250 (1970). 
[47] X. Mao, N. Xu, and T. C. Lubensky, Phys. Rev. Lett. 

104, 085504 (2010). 
[48] C. S. O'Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, 

Phys. Rev. E 68, 011306 (2003). 
[49] L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. Lett. 

95, 098301 (2005). 
[50] C. E. Maloney, Phys. Rev. Lett. 97, 035503 (2006). 



